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Foreword 


Prior to 1990, knowledge of the orbital debris environment could be broken down into two size regimes. For low 
Earth orbiting objects larger than about 10-20 cm diameter, the U.S. Space Command maintains a catalog which 
includes debris as well as intact satellites and rocket bodies. The 10-20 cm diameter limit is the result of the 
wavelengths (~70 cm or longer) and sensitivities of the radars used to detect and track the objects. Knowledge of 
debris smaller than about 0.1 cm comes from analysis of surfaces that have been exposed to space for a period of 
time and then returned to the Earth. The practical size limit of these data result from limitations of the surface area 
exposed and the length of time of the exposure. Models of the orbital debris environment had to interpolate between 
these two size limits (0.1-20 cm). In the late 1980s, NASA needed a better characterization of the orbital debris 
environment to sizes as small as 1 cm diameter in order to verify the design of the Space Station. Under a joint 
agreement with the U.S. Department of Defense (DoD), NASA began using the DoD-funded Haystack radar to fill 
this need. In return, NASA paid for the construction of the Haystack Auxiliary (HAX) radar. 

The Haystack radar is a high-power, X-band (3-cm wavelength), monopulse radar with very high sensitivity. NASA 
began collecting statistical data (detection, but no cataloging) on the orbital debris environment using Haystack in 
August 1990. Since that time it has been the primary source of knowledge of the debris environment in the critical 
size range from 1-20 cm. As such, NASA has striven to continually evaluate and improve the quality of the data and 
the inferences drawn from the data. Several peer review panels have been convened over the years to review the 
Haystack project. The current panel headed by Dr. David K. Barton was convened to answer four specific questions 
dealing with the number and type of observations that have been made and the uncertainty limits which can be placed 
on the resulting size distributions. 

NASA has been converting radar cross section to size based on an empirically derived size estimation model. This 
has been referred to as the “direct” method in the current report. The report has pointed out that, in some size 
ranges, the size distribution is biased. NASA has long realized that this method was biased and, in 1993, began 
using a “weighted” size distribution to correct the bias. One of the panel members. Dr. David Brillinger, has 
developed a more rigorous procedure, presented in Appendix A, which both computes the uncertainty in the size 
distribution and corrects the bias. NASA is in the process of implementing this procedure. 
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1. Summary 


The Haystack Orbital Debris Data Review Panel was established in December 1996 to consider the 
adequacy of the data on orbital debris gathered over the past several years with the Haystack radar, and 
the accuracy of the methods used to estimate the flux vs. size relationship for this debris. The four 
specific issues addressed in the terms of reference [1.1] for the Panel were: 

1 . The number of observations relative to the estimated population of interest 

2. The inherent ambiguity between the measured radar cross section (RCS) and the inferred physical 
size of the object 

3. The inherent aspect angle limitation in viewing each object and its relationship to object geometry 

4. The adequacy of the sample data set to characterize the debris population’s potential geometry. 

Further discussion and interpretation of these issues, and identification of the detailed questions 
contributing to them, are discussed in Section 2 of this report. Section 3 discusses the current status of 
the measurement and analysis program being conducted by the National Aeronautics and Space 
Administration (NASA). Sections 4 and 5 discuss the statistical aspects in use of the Haystack radar 
data, with comments on the issues raised. The conclusions of the Panel with respect to these issues, and 
two recommendations are presented in Section 6 and repeated at the end of this summary. 

A previous study conducted in 1991-1992 by a radar-oriented panel that included two of the present 
panel members had addressed essentially the same issues at the beginning of the Haystack data collection 
effort, in an attempt to estimate how well that radar’ s data would serve as an input to the orbital debris 
flux estimation problem, and to make recommendations aimed at maximizing the usefulness of the 
collection effort. Many of the recommendations of that earlier panel [1.2] were included in the 
procedures adopted by NASA and the Haystack operating personnel in the collection program. In 
particular, the recommendations concerning radar calibration, operating procedures, and selection of data 
to minimize error in measured RCS were addressed in detail, and the resulting error analyses were 
carried out and documented with great care, leading to an estimated accuracy of three to four decibels 
rms in RCS data for a single observation (see Section 3.1 and Appendix A of this report). Recommended 
steps requiring hardware modifications to increase the size of the data sample proved too difficult or 
expensive to undertake. Instead, the sample size was increased by extending the collection period over 
five years and many thousands of hours. 

The remaining issues are concerned primarily with the problem of converting the measured RCS values 
into a model of flux vs. object size (and ultimately weight or other measure of hazard). These issues 
were placed in sharp focus by a draft report [1.3] on a study that personnel of the Los Alamos National 
Laboratory performed for the U.S. Air Force Scientific Advisory Board. It became apparent to the Panel 
that several of the criticisms in this draft were the result of misunderstandings of the actual process used 
in the NASA modeling effort and assumptions underlying that effort. To assist in eliminating these 
misunderstandings, Section 3.2 of this report describes the process by which RCS measurements are 
converted to object size. 
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The remaining issues involve the accuracy of the flux vs. size estimates obtained through use of the 
Haystack data. To address these issues three experts in statistics were appointed to the Panel, and their 
study has led to the following conclusions and recommendations. 

1.1 Issues 

1.1.1 Issue 1: The number of observations relative to the estimated population of interest 


Relative standard error of flux 


Coefficient 
of variation 



Observation time (hr) 

Includes extra Poisson variation multiplier of 1.16 


Figure 1.1. Typical plots of coefficient of variation of the flux (the standard error a of the flux divided 
by the flux) vs. total hours of obseryation for objects in differen t size intervals. To obtain the standard 
error of the flux, multiply the result on the y-axis by the mean flux over the given time interval. 

The number of observations relative to the estimated population is a statistical question and the panel has 
taken a statistical approach to answering the question. An analytical procedure has been developed for 
calculating the number of observations required to estimate the target populations with prestated 
confidence bounds. The model takes into account known sources of error and assumes a family of 
distributions. 

Conditioned on a number of assumptions which are discussed in Appendix A, curves have been produced 
indicating the confidence interval for the population estimate as a function of the total number of 
observations available. Figure 1.1, based on a subset of the Haystack data collected in 1994 and 
representing 840 detections over a 97.22-hour observation period, shows the sort of behavior that might 
be expected as the number of observations increases. The curves in Fig. 1.1, parameterized as a function 
of debris particle size, show that confidence in the population estimate improves with increasing numbers 
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of observations. These curves are illustrative of the type of results that can be obtained, but should not 
be taken as representing the actual values for the entire Haystack data set. 

The details of the confidence interval curves will change as the population upon which they are based 
changes. However, a statistical procedure, explained in Appendix A, is available for NASA to produce 
curves based on other data sets. 

1.1.2 Issue 2 : The inherent ambiguity of the measured dBsm and the inferred physical size 

In the case of the conducting sphere, a given RCS in dBsm can be generated by as many as three different 
sphere diameters, producing ambiguity. The NASA size estimation model (SEM), although it follows the 
curve for average sphere RCS, was derived from many measurements at different aspect angles and 
frequencies on a sample set of 39 objects. The resulting translation of RCS to size is subject to 
uncertainty (not ambiguity), described by a distribution of size values about the mean, for a given RCS. 
Flux vs. size curves can be derived from Haystack data, along with estimates of their uncertainty, 
notwithstanding the variability of RCS with size, shape, and viewing angle. These distributions are taken 
into account in the development of the uncertainty bounds discussed under Issue #1 . 

1.1.3 Issue 3 : The inherent aspect angle limitation in viewing each object and its relationship to 

the object’s geometry 

A Haystack observation of an object is made at unknown aspects, representing samples from a possible 
47t steradians of angle about the axis of the object. The conversion of measured RCS to size, via the 
SEM curve, is based on the average overall aspect angles and the distribution about this average caused 
by object shape and aspect angle variation. The use of this distribution in the derivation of confidence 
limits takes into account the fact that RCS values from only a limited number of aspect angles (typically 
one) are observed in a Haystack measurement. 

1.1.4 Issue 4 : Adequacy of the sample data set to characterize the debris population’s potential 

geometry 

The sample data set derived from a high-velocity impact experiment produced a distribution of shapes 
that is quite similar to the distributions of shapes obtained from other high-velocity impact and explosion 
experiments. Without direct physical sampling of the objects on orbit, this is the most useful evidence 
that the measured data set is representative of objects on orbit. Therefore the measured data set can 
properly be assumed adequate to characterize the space debris population’s potential geometry. 

1.2 Recommendations 

1 . A succinct summary document is needed which describes the debris problem, the Haystack 

measurement program, the procedure for converting RCS to size, and the resulting debris flux results 
obtained to date. This information is embedded in several existing reports, but needs to be integrated 
into a single document for use by those nonexperts in fields of radar and space debris. 
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2. Continue to make use of modern statistical techniques in addressing the interpretation of the 
Haystack data. 

1.3 References for Section 1 

[1.1] Loftus. Jr., J. P., “Data Processing for Haystack Radar Observations,” presented to Haystack 
Orbital Debris Data Review Panel, 10 Dec. 1996. 

[1.2] Barton, D. K., et al., “Orbital Debris Radar Technology Peer Review,” JSC Report, March 1992. 

[1.3] Canavan, G. H., Judd, O’D. P., and Naka, F. R., “Comparison of Space Debris Estimates,” Los 
Alamos National Laboratory Report LA-UR-96-3676, October 1996 (submitted as draft for 
USAL SAB report). 


2. Terms of Reference 

The panel was charged to address the four issues listed in Section 1 . The panel addressed each of these 
four issues from the following perspectives. 


2.1 Issues 

2.1.1 Issue 1: The number of observations relative to the estimated population of interest 

The question posed here is simply whether or not enough measurements were taken to adequately 
characterize the population of concern. The answer to that question is the classic one — more 
measurements would be better! However, given the effort and expense involved in taking the 
measurements and reducing the data, a more detailed answer is essential. 

Lirst, observations made to date consist of radar returns from space debris particles larger than 1 cm at 
altitudes between 350 km and 2000 km trapped in orbit at inclinations above 28° [2.1], The coverage 
across this range is not uniform and can be characterized by the following windows as a function of 
inclination [2.1, Lig. 4.1-2]. 


Beam Orientation 

10° South Staring Data 
20° South Staring Data 
75° East Staring Data 
90° Vertical Staring Data 


Altitude Range 

350 km - 700 km 
350 km - 1200 km 
350 km - 1200 km 
350 km - 2000 km 


Inclinations Observed 

Above 28° 

Above 28° 

Above 42.6° 
Above 42.6° 


Note that the 10° above horizontal, due south, pointing direction of the beam will begin to intercept an 
orbit with a 28° inclination at 530 km altitude. Thus the sampling is far from uniform across the ranges 
identified above. (However, this position and altitude are adequate to reach Space Shuttle orbits.) 
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Given the highly non-uniform sampling of the potential debris environment provided by the Haystack 
geometry and performance parameters, the research team has chosen to use the NASA engineering 
model, ORDEM96 [2.2], to predict the measured environments, and compare the actual measurements 
with this prediction. Since, in the end, a validation of the ORDEM96 model is the puipose of the 
measurements, this is a reasonable and efficient procedure. 

To answer the question of the number of observations relative to the population of interest, the target 
population of interest must be better defined. The primary objective of the Haystack measurement 
campaign is to refine the risk estimate for damage to the International Space Station (ISS) due to space 
debris particles in the range of 1 cm to 10 cm equivalent diameter. ISS will orbit at an altitude of about 
400 km and an inclination of approximately 51° [2.3]. It will be exposed to debris particles of all 
inclinations within this altitude range. Since the majority of estimated debris particles have inclinations 
of greater than 28° [2.2], the measurements made by Haystack do address most of the population of 
interest as it exists today. The confidence that can be placed in how well they measure this population 
depends on the relative accuracy of the reduced measurements. The factors affecting this accuracy are 
addressed in the body of this report. 

2.1.2 Issue 2: The inherent ambiguity of the measured dBsm and the inferred physical size 

The electrical signal the Haystack radar receives when it interrogates a piece of space debris is a function 
of the RCS of the object interrogated. At a given frequency the RCS is a function of the size, shape, 
orientation, and material properties of the object. The determination of the size, shape, and mass of the 
object from a single measurement of its RCS is an impossible task. The estimation of the distribution of 
the sizes of a sample of space debris objects appears to be more tractable. By observing multiple objects 
with probable random orientations and including the physical laws that must govern their aggregate 
behavior, an estimate of the distribution of sizes is significantly more constrained. Objects in this size 
range are likely to have come from breakups. Both explosions and collisions appear to follow power law 
frequency distributions for these sizes. The research team has developed an SEM to convert the 
measured RCS to an inferred equivalent object diameter based on measurements made on debris 
available from hypervelocity impact tests conducted at Arnold Engineering and Development Center 
(AEDC). The adequacy of this process depends on 

1. the ability of the collected debris to represent actual space debris, 

2. the utility of the equivalent diameter concept for characterizing risk, and 

3. the methods used to propagate uncertainty in the estimation process. 

Each of these elements of the estimation process will be commented on in the body of the report. 

2.1.3 Issue 3: The inherent aspect angle limitation in the view of each object and its relationship 

to the object’s geometry 

Any space object detected by the radar beam is in the beam for less than Vi second and the likelihood of 
known multiple detections of the same object is extremely small for the objects of concern. Therefore, 
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the single RCS measurement made during this time will depend on the radar-object orientation at the 
time of measurement assuming modest object rotation rates. This provides an additional statistical 
uncertainty that must be addressed in the size estimation process. It is discussed further in the body of 
this report. 

2.1.4 Issue 4: The adequacy of the sample data set to characterize the debris population’s 
potential geometry 

The sample data set used to build the SEM was composed primarily of objects collected as a result of 
hypervelocity impact tests on satellite mockups. How representative these objects are of actual space 
debris is an open question. Since it is probably not possible to collect a representative sample of space 
debris objects in the size range of interest from space, and measure their RCS, this type of modeling 
inference will have to be used. Analytic tools for, and methods of, characterizing representativeness will 
be discussed in the body of this report. 

2.2 References for Section 2 

[2. 1] Stansbery, E.G. et al., “Haystack Radar - Measurements of the Orbital Debris Environment: 
1990-1994,” JSC-27436, April 1996. 

[2.2] Kessler, D. J. et al., “A Computer-Based Orbital Debris Environment Model for Spacecraft 
Design and Observation in Low Earth Orbit,” NASA TM 104825, November 1996. 

[2.3] Cleghorn, George et al., “Protecting the Space Station From Meteoroids and Orbital Debris,” 
National Research Council, Washington D.C., 1997. 

3. Current Status of NASA Process and Data 

3.1 Gathering of Radar Cross Section Data 

The basic procedures for collecting debris data using the Haystack radar and converting those data to 
calibrated RCS values were evaluated by a previous panel [3.1], and so the details are not repeated here. 
Instead, this paper provides a brief description of the procedure, and then discusses uncertainties in the 
RCS values obtained, as those uncertainties ultimately affect the confidence which can be placed in the 
derived size distributions. 

3.1.1 Calibration 

Haystack utilizes amplitude calibrations and monopulse error voltage calibrations to derive absolute RCS 
values from measurements of power received by the radar. The amplitude calibration of the principal 
polarization (PP) channel relates the power in that channel to the RCS of a calibration target at known 
positions in the radar beam. This is accomplished by tracking a large orbiting sphere (typically having a 
1-m 2 RCS) whose orbital elements and size are well documented. Spiral scans around the sphere are 
used to map the antenna pattern amplitude as a function of beam position and to calibrate monopulse 
angle error voltages in the elevation and traverse channels. In effect, the amplitude calibration procedure 
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provides a calibration constant which relates a particular RCS at a given range and at the peak of the 
radar beam to a certain number of analog-to-digital (A/D) counts at the radar output. That relationship 
between A/D counts and RCS fixes a single point on the radar receiver transfer function. Mapping of the 
complete transfer function is accomplished through an injected signal calibration. In later processing, 
knowledge of the calibration constant, the form of the transfer function, and the negative fourth-power 
relationship between received power and range allows different RCS values than that of the calibration 
sphere at different ranges than the calibration range to be computed. 

Expected RCS values for debris particles of interest lie in the general range from -40 to -20 dBsm. Thus 
the calibration sphere is two to four orders of magnitude larger in RCS than the signals of interest. To 
remedy any concern in the large difference between calibration sphere and target RCS and to evaluate the 
accuracy of the size algorithm as applied to spheres, additional calibrations have been accomplished 
using much smaller spheres in the Orbital Debris Radar Calibration Spheres (ODERACS) experiments. 

In the first experiment, two 5.08-cm-, two 10. 16-cm-, and two 15.24-cm-diameter spheres were released 
and allowed to fly through the Haystack beam exactly as in normal data collection. Reported results 
provided RCS levels within 1 dB of the correct values [3.2], 

A specific calibration concern raised by the earlier Peer Review Panel was addressed in the ODERACS II 
experiment where, in addition to spheres, two 13.348-cm dipoles and one 4.42-cm dipole were released. 
Before that time, there had been no end-to-end calibration of the opposite polarization (OP) channel, as 
the calibration spheres provide only a PP return. Dipoles provide equal return in the PP and OP 
channels, and so a simultaneous measurement of dipole return in both channels, combined with absolute 
calibration of the PP channel, allows calibration transfer to the OP channel. Cress, et al. [3.3], provide an 
analysis of the calibration of the OP channel based on the dipole measurements. The PP-to-OP RCS 
ratios varied between 0.983 and 1.038, indicating that the two channels are well balanced. 

3.1.2 Radar Data Collection 

X-Band radar data are collected with Haystack in a staring mode. Four channels are processed in the 
radar: PP sum, OP sum, PP traverse difference, and PP elevation difference. Data from all four channels 
are coherently converted to a 60-MHz intermediate frequency, filtered to 1-MHz bandpass, further down- 
converted to 5 ± 0.5 MHz, and then digitized at a rate of 20 MHz using a 10-bit digitizer. In-phase (I) 
and quadrature (Q) data are created at a 5-MHz sample rate, and then thinned without averaging to a 
1-MHz rate. Using about a 40% range overlap, the I and Q samples are fast Fourier transformed (FFT) to 
the frequency domain. Complex FFT data for each channel are sent to a memory buffer containing data 
for the previous 12 to 20 pulses. To minimize the archiving of data with no detections, a noncoherent 
12-pulse running sum of the PP sum channel data is maintained, and only when a threshold is exceeded 
are the spectral data for all four channels permanently recorded to tape. The recording threshold is 
intentionally set lower than allowed in subsequent processing to ensure that no usable data are missed. 

3.1.3 Data Processing 

Subsequent data processing accomplished at NASA/JSC converts Haystack-recorded data to detections, 
and from the detections provides RCS estimates. The intent is to have in place an automatic processing 


7 



system, the Orbital Debris Analysis System (ODAS), but manual intervention and operator judgment are 
currently required to ensure only good data are output [3.4], 

The first step in ODAS is an application of a more stringent threshold, again based on the PP sum data. 
For this processing, the signal-to-noise ratio (SNR) calculated in each Doppler bin is adjusted for the 
nonuniform noise floor out of the bandpass filter. Detections again are based on a noncoherent running 
sum of data from 12 sequential pulses, but with a higher threshold than is required for data recording. 
After detection, data from the range cell in which the detection is declared and the two adjacent range 
cells are converted to the time domain, and concatenated to remove the overlap. The resulting data are 
converted back to the frequency domain using a 4096-point FFT and a cross correlation is performed 
against the FFT of a simulated transmit pulse to determine the range to the target. 

The need for the next step in processing is of some concern. Currently, there are phase drift problems in 
Haystack which manifest themselves in monopulse voltage ratios whose phase is neither the 0° or 1 80° 
values which would be expected. To remedy the immediate problem, voltage ratio phases are plotted, 
and then an adjustment is made to the phase, when necessary. MIT/LL is exploring the drift problem, 
and in the presentation to the panel, the data processing team stated that a fix to the problem is imminent. 
It is important that the problem be corrected. A concern of the earlier panel was that contamination of 
the PP difference channel by OP return could bias the apparent angle error, thereby providing inaccurate 
beamshape corrections. If the monopulse voltage ratio phases were routinely 0 or 180°, as expected, 
phases away from those values might be indications of OP channel contamination of angle error data. To 
alleviate potential problems with OP cross-talk signals, the earlier panel recommended using only data 
with relatively large PP/OP ratios. NASA rejected this recommendation because low PP/OP ratio targets 
might be a separate class of targets whose omission might bias distributions of sizes and orbits [3.5]. 
Monopulse voltage ratio phase could be an additional discriminant upon which to base angle error 
validity, if phase drift problems are corrected. 

The final step in the data processing chain leads to average RCS estimates for the debris. In that step, the 
path of the debris particle across the beam is estimated, and adjustments to the measured RCS for each 
data point are made based on the calculated location of the data point within the radar sum beam. Data 
points which meet the criteria for “good” data are averaged to provide the PP and OP RCS values which 
go into the sizing algorithm. 

The path through the beam is estimated using a weighted-least-squares approach. The apparent 
positions — based on spiral scan-derived monopulse calibrations — of the data point with the highest SNR, 
along with the adjacent data points, are calculated. A weighted-least-squares linear path through the 
beam is computed based on the three points, where the weighting is proportional to each point’s SNR. 
Slope and intercept uncertainty terms are also calculated and used to judge the remaining data points. 

The algorithm steps both directions in time away from the three initial points and evaluates each 
subsequent point in terms of whether it falls within the calculated uncertainties. As each point is added, 
the best fit line is recomputed, as are the uncertainties. When a point falls outside the uncertainties, the 
algorithm is terminated in that direction. The panel comments in Section 5 upon possible improved 
procedures to provide better estimates of particle path across the beam. 



Based on the position provided by the best fit line, the amplitudes of the PP and OP sum channels are 
corrected using the map of antenna gain versus position provided by the spiral scans around the orbiting 
calibration sphere. An implicit assumption is made that the PP and OP pattern shapes are identical, as 
the OP pattern cannot be measured using the spiral sphere scans. All points used in computing the best 
fit line which also require less than a 6-dB amplitude modification are included in the average to 
determine the RCS, which is the input to the size estimation routine. 

3.1.4 Error Sources for Individual RCS Measurements 

The data collection and processing procedure has within it inherent error sources which have both 
random and bias components. Major sources of these errors are uncertainties in calibration, propagation 
uncertainties (atmospheric attenuation and refraction effects), signal processing effects (principally 
Doppler straddling losses), errors in the best fit path through the beam, and errors due to additive noise. 

The previous panel, based on experience, estimated a 0.8-dB-bias error and a 2.1-dB root-sum-square 
(rss) random error in radar calibration. A recommendation was made that these values be replaced by a 
more careful calculation using specific Haystack parameters and measurements. Subsequently, Pitts (in 
Stansbery, et al., [3.5]) completed an analysis of measurement error which for the calibration term used 
an MIT/LL estimate of a la random error in the calibration constant of 0.6 dB. In addition to the random 
calibration constant error, bias errors due to atmospheric attenuation were estimated to be -0.7 dB at 10°, 
-0.4 dB at 20°, -0.2 dB at 30°, and -0. 1 dB at 90°. These combined errors provide a comparable case to 
the panel estimate and give a worst-case error which has a 0.6-dB random component and a 0.7-dB bias. 
Details of the analysis were not provided, and so it is not clear if all the error sources (e.g., range 
uncertainty, transmit power fluctuation, pattern propagation factor, etc.) considered in [3.1] were 
included in Pitts’ analysis. 

Errors introduced in signal processing should be principally uncertainties in gate straddling loss. The 
processing step which combines the range bin of the target and the two adjacent range bins should reduce 
range straddling loss to a negligible value. Pitts provides an estimate of straddling loss correction error 
of 0.6 dB for moderate -hit single -pulse SNRs, and this is an unbiased error. 

The largest uncertainties, both random and biased, can arise out of the beamshape correction procedure. 
The previous panel demonstrated that errors in estimation of the particle’s position in the beam lead to 
errors whose magnitudes are dependent on the SNR and the actual position of the particle within the 
beam. Error bias magnitudes and random component standard deviations were shown to increase rapidly 
outside the nominal radar main beamwidth. Based on that finding, a recommendation was made to use 
only data falling within the 3-dB one-way (6-dB two-way) beamwidth of the radar. NASA subsequently 
accepted that recommendation, which was implemented in data processing. Based on the 6-dB two-way 
limits, Pitts estimated the beamshape correction errors for moderate SNR targets to be about a +0.2-dB 
bias error and 2. 3-dB la random error. 
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Pitts also provides his estimates of la bias and random uncertainties due to all causes. As a function of 
elevation angle, they are: 


10° -0.5 dB bias and 2.5 dB random 

20° -0.2 dB bias and 2.5 dB random 

30° 0.0 dB bias and 2.5 dB random 

90° +0. 1 dB bias and 2.5 dB random 


The values given are for moderate (i.e., around 13 dB) SNRs. Note however, that errors will be a strong 
function of the SNR. Table 3.1 provides expected values of single-pulse SNR for a combination of target 
sizes, elevation angles, and ranges. Calculations are based on the radar parameters from [3.2, Table 1] 
and the size estimation model curve which, at 10 GHz, provides a mean RCS value of approximately 
-37 dBsm for a 1-cm effective-sized particle and -21.5 dBsm for a 10-cm particle. The 10° elevation and 
90° elevation cases have been chosen, as they represent the shortest and longest detection ranges that 
might be expected and thus define the limits on expected SNR. Similarly, calculations are shown for the 
particle at the peak of the beam and at the -6-dB two-way point in the beam, as those also represent the 
extremes. Differences in atmospheric attenuation between the two angles have been ignored. 


Table 3.1 Calculated SNRs for Various Elevation Angles and Debris Altitudes 




10-cm debris particles 

1-cm debris particles 

Elevation 

Debris 

SNR (dB) at 

SNR (dB) at 

SNR (dB) at 

SNR (dB) at 

(deg) 

Altitude (km) 

beam peak 

-6 dB point 

beam peak 

-6 dB point 

90 

350 

53.6 

47.6 

38.1 

32.1 

90 

1350 

30.2 

24.2 

14.7 

8.7 

10 

350 

30.8 

24.8 

15.3 

9.3 

10 

700 

22.0 

16.0 

6.5 

0.5 


In all cases for the 10-cm debris, calculated mean SNRs are above 15 dB, and so errors should be no 
worse than Pitt’s analysis. However, for the smaller panicles, that is not the case. Four of the eight cases 
shown have SNRs below 10 dB, and values go as low as 0.5 dB. 

Table 3.2, modified from a table in [3.1], provides expected bias and rss random errors based on SNR 
and actual position in the beam. Terms included are calibration and processing errors, errors due to 
additive receiver noise, and errors due to monopulse angle uncertainty. Pitts’ value of 2.5 dB has been 
used to characterize the calibration and signal processing random component of the error. As Pitts’ 
estimated bias errors are all small, no bias error has been assumed due to calibration and signal 
processing. 
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Table 3.2 Estimated Total Errors After [3.1] 



Beam center 

-6 dB two-way point 

SNR (dB) 

Ibias (dB )l 

la random (dB) 

Ibias (dB)l 

la random (dB) 

25 

0.01 

2.55 

0.01 

2.58 

15 

0.17 

2.95 

0.22 

3.19 

5 

1.97 

6.08 

2.32 

7.22 


For a 25-dB SNR and either beam position, the total errors are dominated by the assumed calibration 
errors. At even moderate SNR values (e.g., 15 dB), the predominant random error source is still 
calibration error, although additive noise errors and inaccurate beamshape corrections due to monopulse 
errors are starting to become important. At low SNRs, the uncertainties are dominated by noise and 
monopulse errors, as calibration and signal processing errors have become a small part of the total error. 
Note the significant bias, as well as the large random error components for low SNR values. 

At the 5-dB SNR level, the la errors due simply to additive noise can provide an apparent SNR of as 
high as 8.9 dB or as low as -2.2 dB. Noise provides both a spread and a significant asymmetry around 
the 5-dB nominal value. This contrasts with the 13-dB case considered by Pitts, where the la uncertainty 
due to additive noise is only 10.8 dB to 14.8 dB. In assessing the confidence intervals on population, the 
change in error bounds with changing SNR should be included. Thus, the panel recommends that a 
careful analysis be conducted of measurement uncertainty as a function of SNR, and that the results be 
folded into the ongoing analysis concerning the confidence which can be placed in the derived particle 
distributions. 

3.2 RCS to Size Conversion 

The current procedure for converting a measured RCS to a debris particle size is relatively straightfor- 
ward. An SEM has been developed based on a collection of fragments that have been measured 
physically and radiatively over the range of relative sizes appropriate for Haystack measurement of 
debris in the 1- to 10-cm size range. Given a RCS measurement, the SEM will assign a unique value to 
the debris object from which it was derived. A detailed discussion of the process follows. 

3.2.1 Effective Diameter 

The size of a particle is defined by its “effective diameter.” The effective diameter is the average of three 
orthogonal thickness measurements. The first measurement is taken along the longest dimension of the 
particle. The second measurement is taken in a direction perpendicular to the first and optimized to 
obtain the largest value possible. The third measurement is made in a direction orthogonal to the first 
two. These three measurements are then averaged to obtain the effective diameter. 

This procedure will of course yield the diameter of a sphere if the debris particle is a sphere. It will also 
permit calculation of the correct volume for a sphere if the object being measured is a sphere. If the 
object being measured is a right circular cylinder, it will overestimate the volume by approximately 38% 
if the height -to-diameter ratio is 1 .0. If the cylinder is stretched into a rod, the approximation to the 
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volume gets better, reaching a ratio of predicted spherical volume to actual volume of 1 .0 about when the 
height-to-diameter ratio hits 2.0 and then diverges to an overestimate of a factor 4 when the height-to- 
diameter ratio reaches 10. As the cylinder is pancaked to a disk, the ratio of spherical volume to actual 
volume hits 1.0 at about a height to diameter ratio of 0.5 and then diverges reaching 2.25 at a height to 
diameter ratio of 0.1. The effective diameter is a useful estimate of size. 

3.2.2 Mapping RCS to Effective Diameter 

Once a procedure has been defined for obtaining the size of an arbitrarily shaped object, it is possible to 
develop a methodology for mapping this size into a RCS, or for performing the inverse and mapping a 
RCS into a size. Figure 3.1 shows the variation in RCS, a normalized to wavelength X, as a function of 
normalized diameter for conducting spheres. The values calculated from Mie scattering theory are 
shown as a solid curve, approximated in the Rayleigh region by a = 9n 5 d 6 /4X 4 . When the sphere 
diameter exceeds 0.2X the RCS oscillates about the optical approximation o = (Tt/4fv 2 , the oscillations 
damped to less than ±1 dB for s > 2X. 



Normalized sphere diameter, sfk 


Figure 3.1 RCS v,v. diameter for conducting sphere. 

A scaling curve for conversion of RCS to size can be obtained by smoothing the oscillations of the 
resonant region and connecting to the straight-line segments in the Rayleigh and optical regions, as 
shown in Fig. 3.2. 

To estimate this function, the NASA research team has chosen a set of 44 objects that were collected 
from two hypervelocity impact shots on satellite mockups conducted at AEDC and some debris-like 
objects collected by the research team [3.6]. The RCSs of these objects were measured over 4 radar 
frequency bands (2.56-3.92 GHz, 4.12-7.986 GHz, 8.15-12.77 GHz, and 12.92-17.54 GHz) with 
eight steps in the band for the lowest frequency and 16 steps in the band for the other three, and with 
approximately 200 source-object orientation angles. This produced a mean RCS for each object, at each 
frequency, and a distribution of RCS for each object depending on source-object orientation. 


12 





Normalized sphere diameter, sfk 
Figure 3.2 SEM curve based on sphere RCS data. 


The measured data are shown in Fig. 3.3 along with the scaling curve shown in Fig. 3.2. It can be seen 
that the measured data on irregular fragments are scattered approximately ±3 dB from the scaling curve, 
except for a series of outliers 3 to 6 dB below the curve. These outliers represent data from a nonmetallic 
printed circuit board that was included in the measured debris. 



SIZE/HAVELENGTH 

Figure 3.3 SEM curve based on measured RCS data [3.5, Figs. A.4.2-1 and A.4.2-2 ]. 

Now if each effective diameter is normalized by the wavelength of the measuring frequency and the 
corresponding RCS is normalized by the wavelength squared, all of the measurements can be plotted on 
the same curve. When these measurements are averaged at each effective diameter-to-wavelength ratio, a 
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single curve can be obtained that is the SEM. This is the curve plotted in Fig. 3.2 with a RCS curve for a 
sphere for comparison [3.8]. Note that the SEM curve goes over to sphere data for both small effective 
diameter-to-wavelength ratios (Rayleigh region) and for large effective diameter-to-wavelength ratios 
(optical region). Only in the resonance region where traveling waves on the surface of the sphere 
reinforce or cancel does the SEM curve miss the sphere’s response. Of course irregularly shaped bodies 
will have much more complicated responses than the sphere’s. The SEM curve becomes the statistical 
calibration curve for the reduction of the measured data. 

3.2.3 Application to Haystack Data 

The direct approach currently being used to interpret Haystack data then proceeds by taking each 
measured RCS and entering Fig. 3.3 on the left, proceeding to the SEM curve, dropping down to the 
bottom axis, and determining the object’s size as estimated by its effective diameter. The observed 
objects are then sorted by size and plotted cumulatively. The uncertainty in the estimate of number of 
particles at a given altitude and for a given size and larger is determined simply by applying Poisson 
statistics to the number counted. Usually 3a error bars are reported. 

The process as it is currently implemented ignores three major sources of uncertainty. 

(a) To obtain an average RCS at a specific frequency for a particular object from the set tested, all of the 
angular measurements made at that frequency are averaged. In the Haystack situation, only one 
source angle geometry is available for each measurement and so there is a distribution of cross 
sections that must be considered for each size of object. This introduces the first source of 
uncertainty. 

(b) Then, to obtain a single curve for the SEM, the average RCS values for the many objects are 
themselves averaged to obtain the average curve. Plots of this data indicate that this process 
averaged data that had a range of approximately 4 to 5 dB (except for the outliers previously 
mentioned). This is the second source of uncertainty. 

(c) The third major source of uncertainty is introduced by taking data down to an SNR that allows a 
large number of false alarms. The count rate vs. size curve is very steep in this region and so it is 
hard to sort out the number of false alarms and correct for them as the size of particles approaches 

1 cm in size. At the 10° above the horizon position and 500-km altitude, a 1-cm particle is no more 
than 5 or 6 dB above the false alarm threshold. If the spread in the SEM curve is 6 dB, one must 
make an accurate accounting of the particles lost because they fell below the cutoff, as well as an 
accurate accounting of below-threshold particles that were counted. This is a difficult problem, but it 
merits attention because the problems with counting statistics are least significant here. 

All three of these sources of uncertainty are at least addressable with improved quantitative analysis. 

3.2.4 Relating the Calibration Data Set to Actual Space Debris 

A major source of uncertainty that is not addressable with better analysis is how well the 44 objects used 
for the SEM database characterize space debris. To address this source of uncertainty, three data sets of 
fragments were analyzed. The first set of data on fragments was obtained from the three European Space 
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Operations Center (ESOC) experiments performed by Fucke [3.9]. These experiments were explosions 
conducted on a scale model of an Ariane upper stage tank. They were carefully performed and over 98% 
of the original mass was recovered after each of the tests. The second set of data was the fragment data 
on the 39 objects recovered from the AEDC tests that have been used as calibration data for the SEM 
curve. The additional objects that have been added to the calibration set were not considered, as it was 
not clear how representative they might be of any fragmentation event in space. The third set of objects 
came from the Satellite Orbital Debris Characterization Impact Test (SOCIT) tests conducted at AEDC. 
These tests were sponsored by the Defense Nuclear Agency and impacted 150-gram projectiles on a 
satellite and several satellite mock-ups to observe the fragmentation patterns likely. The data used here 
were from the TRANSIT satellite test. For this test over 87% of the original mass was recovered [3.10]. 

Unfortunately it was impossible to track all of the measurements in one data file to the complete set of 
data used to generate the SEM curve. Therefore only the original 39 objects were used as representative 
of the calibration set. 

It is very difficult to say that one set of objects represents a large set of objects without at least sampling 
the large set of objects. However, if the set of calibration objects can be used to represent two other sets 
of objects that were obtained by processes similar to what is expected on orbit, then there is reason to 
believe that the calibration set may be representative. Since the measured quantity obtained by Haystack 
is the radar cross section of the fragments, and the parameter used to characterize an object is its effective 
diameter, the most logical parameter to use to compare the representativeness of one sample to another is 
the ratio of an area presented by an object to the square of its effective diameter. The projected area of 
an object can be estimated as being proportional to the product of any two orthogonal dimensions. The 
effective diameter is calculated as the average of the three orthogonal dimensions. One measure then for 
comparing two sets of objects is to calculate the average of the ratio of the product of the two longest 
dimensions to the effective diameter squared. A second measure is to use the product of the two smallest 
dimensions, and a third measure is to use the product of the largest and the smallest. This has been done 
for each of the data sets and the results are presented in Table 3.3. 


Table 3.3 Comparison of Parameters for ESOC, RCS, and SOCIT Data Sets 



R1 

R2 

R3 

R4 

ln(R1) 

ln(R2) 

ln(R3) 

ln(R4) 

ESOC 

Average 

1.528 

0.365 

0.622 

783.785 

0.404 

-1.158 

-0.577 

6.616 


Stnd Dev 

0.290 

0.183 

0.240 

270.591 

0.203 

0.598 

0.513 

0.300 











RCS 

Average 

1.518 

0.437 

0.667 

243.823 

0.400 

-0.959 

-0.501 

5.326 

Calib 

Stnd Dev 

0.276 

0.202 

0.253 

140.410 

0.195 

0.578 

0.517 

0.625 











SOCIT 

Average 

1.411 

0.448 

0.701 

346.850 

0.332 

-1.048 

-0.464 

4.888 


Stnd Dev 

0.231 

0.266 

0.290 

655.736 

0.165 

0.799 

0.521 

1.351 
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In Table 3.3, the quantity R1 represents the ratio of the largest dimension times the next largest divided 
by the effective diameter squared. R2 represents the same ratio using the two smallest dimensions, and 
R3 represents the largest times the smallest divided by the effective diameter squared. R4 is the ratio of 
the effective diameter squared divided by the mass of the particle in units of mm 2 /gram. The last four 
columns are simply the averages and standard deviations of the logarithms of the ratios. 

It is apparent that the average dimensional ratios for all three data sets are very close to each other. 
Therefore, if the materials are very similar, the average measured RCSs ought to be very similar. For this 
problem, the RCS calibration data set is representative of the other two and probably the material in 
space. Therefore, for obtaining effective diameters for objects observed with Haystack, the RCS 
calibration set should be representative. 

However, the fourth ratio is dramatically different between the three data sets. Only a very small part of 
this difference is due to the possibility that this ratio will vary slightly with size of object. The major 
difference is due to the different materials in the three data sets. The ESOC data is mostly aluminum. 

The RCS calibration data has a lot of stainless steel in it, and the SOCIT data set has a lot of phenolic 
material in it. The composition of the actual material in space is an unresolved question. Thus the RCS 
data set is probably representative of the shape of objects in space, but a great deal more analysis is 
required to determine how representative it is of the mass of debris pieces in space. It is really the mass 
of an object that determines its penetrating power, and the conversion from effective diameter to mass 
will need to be addressed in shielding studies. 

3.3 Comparing Measured Data With ORDEM96 

It is very difficult to go from data measured crossing the beam of the Haystack observation volume to an 
equivalent flux of particles orbiting about the earth. And in the end, the purpose of the measurements 
with Haystack is to validate the engineering model for estimating space debris fluxes with the ORDEM96 
model. Therefore the most straightforward approach is to calculate the ORDEM96 particles passing 
through the Haystack beam and compare the resulting count rates with those observed. An example of 
this comparison is Fig. 6 of NASA TM 104825 [3.7], reproduced here as Fig. 3.4. To understand this 
type of comparison and how it validates ORDEM96, a discussion of ORDEM96 is in order. 

ORDEM96 is the NASA engineering model for predicting the debris environment in near earth space. It 
tracks debris by breaking up the ranges of the independent parameters into sub-ranges required to 
describe a debris particle’s trajectory over the earth and calculates an average value for each of the sub- 
ranges. Basically there are six independent variables required to track a particle in inertial space. These 
could be the three components of its radius vector and the three components of its velocity vector at a 
specific time, or they could be the orbital elements (semi-major axis, eccentricity, inclination, argument 
of perigee, right ascension, true anomaly). ORDEM96 assumes that debris particles are scattered 
randomly about their orbits and assumes that any given value of mean anomaly is equally likely. It also 
assumes that right ascension (the point that the particle crosses the equator as it passes from the southern 
to the northern hemisphere) and the argument of perigee (the location of perigee in the plane of the orbit) 
are randomly distributed. All three of these assumptions are reasonable with two exceptions. 
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Particle diameter [cm] 


Figure 3.4 Comparison of ORDEM model for cumulative debris flux with Haystack data [3.7], 

Molniya orbits used by Soviet communications satellites are located with an inclination of approximately 
63° and an argument of perigee in the southern hemisphere. This orbit was chosen because the argument 
of perigee does not precess or randomize with these parameters. Very little is known about the debris 
generated by the Molniya satellites and therefore they are ignored as a debris source. These orbits tend 
to have a very large eccentricity and therefore spend very little time crossing low earth orbits so their 
neglect may not be a major discrepancy. 

Assuming a random distribution of right ascensions for orbits with inclinations slightly over 90° (sun- 
synchronous) may not be reasonable. Sun-synchronous satellites are designed to pass over the same spot 
at the same time of day every day, so radar observations may be correlated with time of day for these 
orbits. This probably is not a problem in using ORDEM96 to predict fluxes on spacecraft, but could 
cause measurements to bunch at certain times. 

This leaves the final three parameters, semi-major axis, eccentricity, and inclination. Since most 
launches to date have been made in the northern hemisphere and it is more efficient to launch in an 
easterly direction to obtain the advantage of the earth’s rotation, most satellites are in orbits with 
inclinations less than 90°. In fact different functional programs have tended to put debris into orbits with 
similar inclinations. Therefore ORDEM96 breaks the inclination variable up into six ranges. They are: 
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ORDEM96 Inclination Boundaries 


0° to 19° 

19° to 36° 

36° to 61° 

61° to 73° 

73° to 91° 

91° to 180° 

Each of these ranges can be correlated with launch sites and particular missions, so that is where most of 
the debris is likely to be. The one exception is the range from 91° to 180°. It is very likely that the range 
from 91° to at most 1 10° is populated, but the range from 1 10° to 180° is very unlikely to contain a 
significant amount of debris due to the difficulty in getting to this range, and the lack of a utility bonus 
for so doing. 

The remaining two parameters are handled by assuming that orbits are either circular or have an 
eccentricity equal to 0.2 with an apogee of 20,000 km. This is a reasonable assumption and adequately 
predicts the fluxes observed on the trailing surfaces of the LDEF satellite [3.7], For circular orbits, the 
semi-major axis is the radius from the center of the earth, and we can derive the semi-major axis for 
elliptical orbits with a fixed apogee and a defined eccentricity. The radius of circular orbits is treated as 
a continuous variable in ORDEM96 so ranges can be defined at the time of comparison with data. 

ORDEM96 also breaks all debris up into components based on its size. The following components are 
defined. 


Component 

Intact objects 
Large fragments 
Small fragments 
NaK particles 
Paint flakes 

Aluminum oxide particles 


Size range 

50 cm <D 

1 cm < D < 50 cm 

200 microns < D < 1 cm 

200 microns < D < 1 cm 

20 microns < D < 200 microns 

D < 20 microns 


This is obviously not a breakdown strictly by size class because not all types of particles exist at all 
inclinations or altitudes. The source of NaK particles is believed to be leaked coolant from a spent 
Soviet RORSAT; it should only exist at one inclination (65°). The aluminum oxide particles appear to 
come from solid rocket burns and should only exist where solid rockets have been used. There are no 
circular debris objects in the less than 19° inclination band, as the debris has primarily resulted from 
satellites launched to geosynchronous orbits. 
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Haystack observes the intact objects and the large fragments above 28° inclination. Intact objects and the 
high end of the large fragments are well measured by the U.S. Space Command Catalog. Its main utility 
is its ability to measure the low end of the large fragment population. 

A comparison can be made with the ORDEM96 model by calculating the number of debris particles 
“colliding” with the Haystack beam during its hours of operation at a given azimuth and elevation. 
ORDEM96 predicts the size of the particles, their altitudes, and their inclinations. The Haystack radar 
measurements give the size and inclination distributions within an altitude band. These two estimates of 
the same population can be compared to quantify how well ORDEM96 predicts the measured data. 
Unfortunately the uncertainties in the measured data tend to be dominated by the low count rates 
observed and the data have to be aggregated in fairly large bins to obtain useful comparisons. 
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4. Ongoing NASA Work 


A basic problem that NASA is studying is the estimation of the debris size distribution given measured 
RCS values. See Stansbery et al [4.1] and Matney [4.2] for descriptions of the problem: Use is made of 
collections of size and cross section data available from laboratory experiments and of recorded Haystack 
cross section observations. In an important step forward, Matney sets up the problem as one of solving 
an integral equation expressing the Haystack measured cross section distribution as a functional of the 
desired debris size distribution and a calibration distribution (see Eq. (4.1) below). Previously an explicit 
functional, as opposed to distributional, relationship has been used to relate cross section to size (see 
Fig. 3.3 above and Section 5.3.3c). The estimate constructed is referred to as the direct estimate below. 
The computations involved (SEMs) are part of NASA’s program SIZEST. 

To begin, some connection between observed RCSs and physical size must be developed from the 
available laboratory data. Connections may take several forms: 

i) . Interpolated: A simple interpolation curve is set up between theoretical curves for the optical and Rayleigh 

regions. 

ii) . Functional: A curve, s = ffz), is obtained from (dimensionless) size and reflectivity values, s and z. 

iii) . Distributional: One estimates the conditional distribution p(zls), of the z values, given an s value. 

Assuming p(z\s) is available, Matney [4.2] sets down the relationship 

m(z) = j n(s) p(zls) ds (4.1) 

connecting a desired distribution n(-) of sizes of objects and a measured distribution m(-) of RCSs. An 
iterative method of solution is proposed for (4.1) and some examples are presented. Practical difficulties 
are mentioned, such as unequal probabilities of detection. The procedure is illustrated for the 
unweighted and the weighted SEM (i.e. both bias-corrected and not). 

The equation (4.1) and its iterative solution are discussed in Vardi and Lee [4.3]. The equation is 
Fredholm of the first kind and is so-called ill-posed (see Allison [4.4]). The ill-posedness leads to 
difficulties in the solution (see O’Sullivan [4.5] and below). Another difficulty with the approach is that 
since no error distribution is involved, there are concerns as to the efficiency of the estimates. 

An error analysis of the solution is essential for the flux estimates. Matney [4.2] also presents some 
ideas for obtaining estimates of variance. He proposes to develop statistical error bounds for the solution 
via the method of propagation of error. This procedure is based on linear approximations to the functions 
involved and requires the availability of variance estimates of basic quantities. Covariances are assumed 
to be zero. 

Difficulties with an approach based simply on (4.1) include that it is not clear how to incorporate 
supplementary information beyond simply studying subsets of the data, nor is it clear how to stabilize the 
estimates. In particular one can think of other variates to include, such as polarization, the complex (or 
matrix) valued RCS, solar activity, secular, and abrupt temporal variation. 
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A stochastic modeling approach in many cases has the pertinent flexibility and that is what will be 
focused on below. 

4.1 References for Section 4 

[4. 1] Stansbery, E. G., et al., “Haystack Radar Measurements of the Orbital Debris Environment: 
1990-1994,” JSC-27436, April 1996. 

[4.2] Matney, M. J., “A Method for the Computation and Error Estimation of a Debris Size 
Distribution from a Measured Cross Section Distribution,” JSC Report, 1996. 

[4.3] Vardi, Y., and Lee, D., “From Image Deblurring to Optimal Investments: Maximum Likelihood 
Solutions for Positive Linear Inverse Problems,” J. Roy. Statist. Soc. B55, 1993, pp. 569-612. 

[4.4] Allison, H., “Inverse Unstable Problems and Some of Their Applications,” Mathematical 
Scientist 4. 1979, pp. 9-30. 

[4.5] O’Sullivan, F., “A Statistical Perspective on Ill-Posed Inverse Problems,” Statistical Science 1, 
1986, pp. 502-518. 


5. Detailed Remarks on NASA Work 

This section is written from a statistical standpoint, in contrast with Sections 2 and 3 that have an 
engineering origin. Necessarily the material contains some review of the existing approaches. The basic 
suggestion is that NASA look into and implement a variety of pertinent statistical methods. The practice 
of some of these methods is now fairly well-developed. 

5.1 The Problem 

The problem at hand has been variously described as: 

“To characterize the debris environment in terms of size, altitude, and inclination.” 

“To predict flux (as function of size, velocity, direction, ...) at specified location in low earth orbit.” 

“To define the orbital debris environment.” 

“To statistically characterize the orbital debris population.” 

“To examine how the orbital debris environment changes with time.” 

“To make environmental predictions - flux levels, directionality...” 

“To provide flux onto a spacecraft surface, the collision rate.” 

And, finally, the specific goal “to assign realistic error bars to the flux curve derived from the Haystack 
data.” 

The problem will be taken as the description of the existing debris environment and the prediction of flux 
values with indications of the uncertainty of estimates constructed and the assessment of the validity of 
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the models employed. The material presented has in mind three principal problems: developing a 
calibration procedure relating physical sizes to return radar reflections, employing the data collected to 
get a flux-size curve, and, finally, validating procedures like ORDEM96. 

Other interesting problems include extracting debris families, for example by altitude versus range rate 
plots and the related problem of assessing how many objects populate a given orbit. 


5.2 Recognizable Component Pieces 

NASA has established a structure for their overall work. They have developed associated procedures and 
computer programs. To set the scene, we provide a list of some of the components, via NASA’s names 
for the computer routines, and then in the following sections we provide a description of pertinent 
statistical methods with details laying out the connection with the NASA work. 


NASA Routines: 

detects passages of objects through the radar beam on the 
basis of returned signals. 

fit a line of passage through the radar beam and convert the 
result into orbital elements and cross sections. 

relates the flux of objects passing through the beam to the 
number detected. 

the size estimation model. 

the model used to predict expected flux level and velocities as 
a function of various quantities and relationships; prediction 
may be for a telescope’s field of view or for a spacecraft. 

5.3 Statistical Techniques Pertinent to the Component Problems 

The field of statistics has something to offer the various component areas of the structure referred to 
above. The basic problems and associated statistical methods are next addressed. 

A common element of contemporary statistical analysis is a focus on stochastic models. These models 
have important advantages, including efficiency of analysis, uncertainty computation, assessment of 
goodness of fit, and combination of data sets. Setting a specific model down makes clear the working 
assumptions. There are stochastic model descriptions associated with each of the methods highlighted 
below and there is also an attempt below to note the models latent in NASA’s work. 

The debris environment is both spatial and temporal. Having fixed on a particular range of altitudes and 
inclinations, the probability density of the sizes of the objects may be assumed stationary for restricted 
time intervals and this density is one of the things to be estimated. There is also the temporal aspect 
relating to the rate at which objects are passing through the field of observation. This may be described 
as a count of objects through part of the beam or converted to a flux. Either of these quantities may be 


HSTRIP: 

HCOR, SIZEST and 
ELSET: 

RPM ( radar 
performance model ): 

SEM: 

ORDEM96: 
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multiplied against the size density estimate to obtain rates for objects of some particular size. The rate of 
passage will be changing in response to external quantities (covariates) such as solar activity. 


Section 5.3.1 refers to the fundamental problem of relating RCS to physical size, that is, the problem of 
developing a statistical calibration model. This statistical calibration activity represents separate work 
from that of flux estimation. Its results apply to land-based objects being illuminated by radar as well. 
This problem can be worked on quite separately and one expects the solution to steadily improve with 
time as the results of more laboratory experiments become available. Section 5.3.2 provides some 
comments related to the determination of orbital parameters after a detection of an object in the beam. 
Section 5.3.3 addresses the problem of estimating the size distribution of the objects in space. Section 
5.3.4 refers to the fundamental problem of flux estimation and 5.3.5 and 5.3.6 to techniques for handling 
the difficulties of biased observation and errors in measurement. The problem of uncertainty estimation 
is addressed in Section 5.3.7. 

5.3.1 Statistical Calibration 

To begin, to avoid possible confusion by some readers, it is remarked that the problem of statistical 
calibration is distinct from that of radar calibration. In the statistical case the object is to find a 
functional description so that observations of one form may be effectively replaced by those of another. 
There may be a calibration curve, such as s - g(z), or a calibration distribution, p(s,z). There are 
parametric and nonparametric procedures (see Clark [5.1, 5.2]). Osborne [5.3] provides a review. 

When a satellite passes through the radar beam, it may be detected and a number of quantities measured 
and satellite characteristics inferred. In particular the RCS and orbital parameters may be mentioned. 
(Remember, for later refinements, that these values are themselves subject to measurement error.) The 
time of the observation will be known. 

In the present calibration context, the basic problem is to relate the return signal (RCS) to the size of an 
object detected. One has to be concerned with whether or not the problem actually has a solution. In that 
connection, Levanon [5.4] states that “most targets exhibit different a (RCSs) at different aspect angles 
and at different frequencies,” but goes on to present: 


Table 5.1 Typical RCS Values 


Target 

o (m2) 

Target 

a (m2) 

Insect 

10 5 

Small aircraft 

2 

Bird 

0.01 

Large fighter 

10 

Small missile 

0.1 

Large airliner 

40 

Man 

1 

Car 

100 


These values fit in with common notions of size (although the “small aircraft” versus “car” values give 
pause for some thought). 
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Stansbery et al. [5.5] and Matney [5.6] describe calibration procedures, motivated by their Figures 3.2-1 
and 1-3 respectively (see Fig. 3.3 of this report), which relate RCS and physical size (mean of maximum, 
maximum perpendicular, and final lengths) on the basis of measurements coming from a laboratory 
experiment. Work is with dimensionless size values, s — S / X, and RCS values, z = CJ / X 2 . The 
statistical discussion of this section applies to other measures of size, e.g. mass, equally. 

As mentioned in Section 4, there are three approaches to the problem: simple interpolation between the 
Rayleigh and optical curves, fitting a functional relationship, and working with the joint distribution of 
size and RCS values. The first two approaches are effectively the same and correspond with a statistical 
model 


z = h(s) + £ (5.1) 

In the first case the form of h is simply assumed, in the latter it is estimated from laboratory data, e.g. by 
taking the average of all the observed z values for s in a chosen interval. To proceed next an estimate of 
size, for given z, is derived via 

s = h~\z) = g(z) (5.2) 

The estimated functions h and g are subject to uncertainty being based on a sample of values. By 
repeated laboratory experiments, this sampling uncertainty may be reduced steadily. (It is assumed away 
in the direct approach.) This is not the case for the uncertainty due to £ in (5. 1). This latter quantity has 
two basic variance components, one due to variability between objects and a second due to variation of 
viewing angle. The variation due to viewing angle is ever present and presents a major problem for the 
estimation of the size distribution, as will be detailed in Section 5.3.3. 


In the distributional approach one estimates, for example, the density function p(z\s) of z values for a 
given s value. The estimate may be of parametric functional form. For example Bohannon et al. [5.7] fit 
mixtures of gamma and log-normal distributions to the basic laboratory observations. Alternately it may 
be a nonparametric, e.g. kernel, estimate. Actually, the model (5.1), if considered in full, corresponds to 
a distributional approach with 


P(z | s) = 



[z - h(s)f 
2<5 ' 


(5.3) 


if one is willing to assume that £ has a normal distribution. 

For both the functional and distributional approaches well-established statistical procedures are available 
for obtaining uncertainty estimates of parameter estimates and predicted values. Section 5.3.3 argues that 
the distributional approach is the more appropriate to employ in estimating the size distribution.. 

To end this discussion, it is appropriate to remember that there are the problems of just what is an 
appropriate measure of physical size and also to remind the reader of the cause for concern involved in 
the use of relationships derived in the laboratory for the orbital case. The basis is an intellectual 
argument. For example, astronomers often have to make use of assumptions that are not directly 
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verifiable. They look for consistency of various inferences to gain confidence in the results. If a variety 
of laboratory studies appear to be consistent, then one may gain some confidence in extrapolating to the 
orbital case. In the present case there may also be ways to infer the relative proportions of various 
materials in orbit and make use of this information in the extrapolation. 

5.3.2 Estimation of a Satellite’s Path Through the Beam 

Before addressing the problem of estimating the size distribution of the objects in orbit, basic 
measurements of size-related values are needed. One can detect the presence of individual objects 
passing through the beam. Estimates of the parameters of their paths, particularly altitude and 
inclination, are required for later analyses. Presently NASA applies a weighted-least-squares procedure 
to the data values surrounding a detection employing subjective decision-making to accept or reject in 
particular cases. 

It is worth noting that, supposing that a detection has been made, the tools of robust regression analysis 
can be employed to fit a (straight line) path through estimates of position. Robust procedures are meant 
to automatically and objectively reduce the influence of extreme observations (outliers). Weights to 
handle varying levels of statistical uncertainty are incorporated simply and uncertainty estimates, such as 
standard errors, are part of the usual output. Books by Huber [5.8] and Rousseeuw and Leroy [5.9] cover 
the topic. Anderson-Sprecher [5.10] is a reference to the robust estimation of location. Netanyahu et al. 
[5. 1 1] applies such procedures to detect road segments in images. 

In a variety of analyses the detections will be cross-classified, for example, by altitude with the primary 
concern as expressed in ORDEM96 for objects passing near the intended location of the future space 
station. 

Remember when using the altitudes and inclinations that they are estimates themselves, subject to 
measurement error. 

5.3.3 The Estimation of the Size Distribution 

The concern now is that of obtaining an estimate of the distribution of the object sizes in orbit on the 
basis of RCS values obtained via Haystack. There are in fact a number of approaches to this problem. It 
will be argued that the direct estimate, currently being employed by NASA, is biased (see Section b 
below). However, alternate approaches are available for which the bias will fall to zero for large sample 
sizes. In fact, likelihood-based procedures may be developed and these may be anticipated to be 
efficient. Also, regularized procedures are available to incorporate smoothness assumptions and to 
stabilize the estimation procedure. 

Assume for the moment that statistically the collection of objects being considered is homogeneous, not 
varying substantially through solar activity or as a result of explosions instantaneously creating new 
objects. Thus, remember that the estimated quantities are appropriate but for the particular time period 
over which the data values have been collected or for other periods with similar levels of solar and 
launch activity. 
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However estimated, the distributions of particle sizes, and analytic approximations to them, are crucial 
input to programs such as ORDEM96 and for risk estimation. 

a. Formulation of the problem. 

The problem may be described as estimating for some set of orbits the number or flux of objects with 
size greater than a given size or with size in a particular interval of size values. The problem seems 
usefully split into two parts: a first of determining the overall probability density of the sizes and a 
second of determining either the number of objects with particular characteristics or the corresponding 
flux levels. These last may be handled as multiples of the cumulative distribution function or of the 
density estimates. One can anticipate that, for large enough sample sizes, these last estimates will be 
close to their population values. Counts will remain statistically variable, while rates will stabilize as 
additional data of the same sort are included. 

b. The direct estimate. 

The current procedure for a particular time period and telescope field of view (with restriction to some 
time period and some subset of orbital parameter values) is: Given Haystack observations, 

Z„,n = 1,..., (V , use the relation (5.2) to obtain size estimates s n and then estimate the cumulative size 

distribution via 


1 

N 


I 


H(s n -s 0 ) 


(5.4) 


where Hi ) is the Heaviside function. (His) = 1 for 5 positive and = 0 for s negative.) Unfortunately this 
estimate is generally biased by the nonlinearity of the procedure, biased in a way that collecting more 
data values will not alleviate. 


To see this suppose, as an approximation, that in some interval the function h(s) of (5.1) may be 
approximated by his) = a + /3s , then 


(z„ -&)/$ = s„ + 


s n + — 

" p 


[(q-a) + (P-P)s„ +e„] 

P 


(5.5) 


At the last step it has been assumed that the variability of 0C,(3 is negligible (as will be the case with 
sufficiently many laboratory experiments). From (5.5) it is apparent that the distribution of an estimated 
size s is the convolution of the desired distribution of 5 and of an error term £ / 1) . Acting as if the £ n 
have the same density function f E (•) , 


E{H(s-s 0 )} = iG 


P. 


/ e (£)r/8 


instead of the desired complimentary c.d.f. Probjv > ,s 0 } = G(s 0 ) . For the expected value to be close to 
G(Sq ) , the error £ will need to be small or (3 large. 
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The magnitude of the bias may be investigated by Monte Carlo methods on available data and 
theoretically as follows. Consider the following estimate of the proportion of the objects with sizes in 
the interval (a,b): 


-^E 


(5.6) 


It has expected value 


P. 


\ds 


f e (£)dz 


Assuming that £ is normal with mean 0. variance O” and that O / (5 is not too large, this expected value 
is approximately 



\f\b)-f\a)]<5 2 

2(3 2 


i.e., the estimate is biased and the bias may be estimated from the available data. Examination of the 
figures in Matney [5.6] suggests that the direct estimate of the cumulative count is biased downwards 
The assumption of approximate linearity in z = h(s ) = a + fis might be replaced by one of log z is 
approximately linear in log s . 


There are methods to debias estimates like (5.4) and (5.6) (see Gaffey [5.10] and Zhang [5.1 1]), but their 
specific applicability and effectiveness needs to be investigated for the present situation. 

In summary, the estimate currently being employed is biased and alternate approaches need to be studied. 

c. A grouped estimate. 

An alternate approach is available for determining statistically reasonable estimates for the present 
situation. It is based on the relationship 

p(z ) = I p(zls) f(s)ds (5.7) 

connecting the desired size density /(•) with the directly estimable p(-) and p( T). Specifically, the 
problem may be simply approximated by a Poisson regression analysis and these analyses have been 
often studied. 

Suppose that the observed RCS values have been grouped into intervals I k with the k-th centered at u k 
and containing m k measured RCS values. Suppose the desired range of s values has been broken into 
intervals 7/. Now the distribution of the m k may be approximated by a multinomial with cell 
probabilities 


Pk = L Pmfi 

I 
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setting down a discrete approximation to (5.7). In particular, one could take 


ft =|, f(s)ds 

J Ji 

and 


Pk\i = P(P k \Si) 

for s t the central value of the interval J t . 

Assuming that the p kll are not subject to substantial variability (remember they have been estimated as in 
Section 5.3.2 above), the f, may be estimated by standard generalized linear model techniques 
(specifically Poisson regression) and associated measures of uncertainty are available. These methods 
are described in McCullagh and Nelder [5.14], for example, and various computer algorithms are 
available. Actually the iterative procedure of Vardi and Yee [5.15] is one way of determining the 
maximum likelihood estimate, but the iteratively reweighted-least-squares approach of generalized linear 
modeling appears more flexible. 

In this multinomial case the log likelihood to be maximized may be taken to be that of the Poisson, 

Ik log(A^P„/ ; ) - NYjPkJA (5-8) 

k l I 

where N = '^ m k . It is to be maximized as a function of the /) . 

Once estimates of the /, are available, we may derive estimates of the cumulative distributions by 
summation. In making this procedure operational, the conditional density p(z\s) might be taken to be 
normal, as represented by (5.3) or the functional forms in Bohannon et al. [5.7]. 

One must address the problem of the choice of I k and /, : If too many intervals J , are employed, the 
estimates become unstable. 

This approach has the flexibility of allowing covariates, such as solar activity, launches, and trend to be 
included. For example, one might employ the expression 

(!>«/, ) ex p{A,J 

i 

to correspond to the observation z n € I k that has covariate value x n . 
d. Fixed likelihood-based procedures. 

With p(z) given by (5.7), and 0 denoting a parameter describing the desired size distribution, the 
likelihood based on the z. n data values may be written 

L(Q) = n p(z n ) = Wp(z„ \s)f(s\Q)ds 

i i 
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and 0 then estimated by maximization. For example, /(.s]0) might be taken to be a mixture of power- 
law (or Pareto) distributions as mentioned in Stansbery et al. [5.16]. 

A related estimation approach is provided by nonparametric procedures. These provide a highly flexible 
method to estimate a distribution. One assumes, for example, that 

logf(s\Q) = Q l B l (s)+...+Q I B I (s)-C(e) 

with the B x (•), cubic splines and the dimension I estimated by minimizing the Akaike or 
Bayesian information criterion and then the 0’s by maximizing L(0). Stone [5.17] is a reference to a 
simpler situation. 

e. Regularized estimates. 

An estimate determined in the grouped approach of Section c above is of histogram form involving 
intervals of constancy and interspersed with jumps. Further, the result obtained depends on the choice of 
cells and how the integral of (5.7) is discretized. 

Regularization is a flexible method for determining smooth estimates. It has the additional property of 
making the estimates stable. The estimates obtained by the Poisson regression procedure described 
above may well be unstable if the number of cells for the s variable is taken to be too large. One method 
of regularization is to add a penalty term to the log likelihood. For example, a term 

*•£(/, -/,+,) 2 

might be added to (5.8). The quantity X is called the smoothing parameter. It may need to be estimated 
as well. Cross-validation is one procedure for obtaining an estimate of X . 

The computations of this optimization may again be carried out after a discretization, but they may 
sometimes be conveniently carried out via a variant of the EM (expectation-maximization) method, see 
Green [5.18], He illustrates the procedure specifically for the case of Poisson regression. 

The study of such penalized procedures is related to the general study of inverse problems. This is an 
important subfield of contemporary applied mathematics and statistics (see O’Sullivan [5.19] and 
Tarantola [5.20]). Such problems are particularly difficult when they are ill-posed, that is, the 
corresponding inverse operator is unbounded as is the case for (5.7) (see Allison [5.21]). 

To set up the situation in inverse problem form one can proceed as follows. Introduce a parametrization, 
such as 


f(s) = exp{0(^)} / 1 exp {Q(s)}ds 

for the density. This ensures that /(.) > 0, j f{s)ds = 1 . Following Gu [5.22], (see also Cox and 
O’Sullivan [5.23]), one might now seek 0(-) to maximize 

E log p(z n ) - N log I exp [0( s) }ds + X j 0(s) 2 ds 
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Again A, may either be fixed or estimated. 

Iterative solutions for related problems and procedures for estimating X by cross-validation may be 
found in Wahba [5.24, 5.25], O’Sullivan [5.19], and Gu [5.22, 5.26], The books by Silverman [5.27], 
Tarantola [5.20], Hastie and Tibshirani [5.28], and Wahba [5.29] present various related ideas. 

5.3.4 Flux or Intensity Description and Estimation 

Other problems include extrapolation from the observations of the field of view to other regions of the 
sky and from restricted time intervals of observation to others of interest. Graphs reported variously 
provide cumulative number, cumulative detection rate (number/hour) and flux (number per square meter- 
year) versus size for some population of objects of interest. One of the density estimates of the previous 
section may simply be “multiplied up” to provide such. In the first two cases one uses the overall count 
of objects in the sample. Remember that, for a given altitude for example, the count depends directly on 
the width of the altitude bin employed. For the flux estimate, in the case of vertical staring, one can 
make use of the simple relation 

AA = 2nhAh tan 0/2 

for the cross-section area of the beam at height /?, width Ah and beamwidth 6 . If N objects are observed 
in a time interval of length T, then the overall rate of passage through the beam my be estimated by N/T 
detections per unit time and the overall flux as N / TAA . Alternatively the Radar Performance Model 
may be employed in some fashion. 

The satellite population characteristics are dynamic: There are births and deaths and objects are 
changing their parameter values. Stochastic point and particle processes are basic probabilistic concepts 
used in describing locations and motions of objects. They are described by (conditional) intensity 
functions involving parameter values and explanatories. To the extent that they provide a desired flux 
rate, their description may be seen as the goal of the work. The power of the stochastic point process 
approach is substantial. One reference is Fishman and Snyder [5.30]. 

We now need a process describing the times at which particles with particular characteristics, e.g. size, 
pass through a given field of view. The stochastic description would allow inclusion of covariates such 
as trend or solar activity. Estimates of the parameters of the process may be derived from stretches of 
data. A variety of special stochastic process models are available. 

ORDEM96 may be seen as a procedure whose goal is to predict the intensity function of a spatial- 
temporal marked point process. Its output is directly pertinent to the problem of risk estimation. 
Unknown quantities, e.g. counts of particles with particular characteristics, involved in the problem are 
estimated (see the Appendix of Kessler et al. [5.3 1]). The time interval dt of ORDEM96 is a year and the 
program further makes assumptions concerning the future behavior of some quantities, e.g. solar activity, 
rate of production of new debris, and launch rate. 
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5.3.4 Observational Biases 


There is the added difficulty of unequal probabilities of observation. Stansbery et al. [5.16] apply 
weights to an estimated calibration curve and a modified curve; their Figure 5.0-4 is used in later 
computations. Matney [5.6] lists three sources of bias: false alarms, measurements near the lower limits 
of possible observation, and objects not passing through the center of the beam. The RPM evaluates the 
probability of detection and the edge-to-edge distribution necessary for correcting the data. 


In general terms one can proceed as follows. Let the “true” distribution of the z values be p(z). Let the 
probability that a satellite with RCS = z is detected be Jt (z). This probability will also depend on 
inclination for example, but that will be ignored for the moment. The function n(z) will be near 1 for 
“most” z values and decrease with z. The density of the available RCS values is now 


m(z) = 


n(z)p(z) 

J 7 i(z)p(z)dz 


The likelihood of the data, zi , ■■■ , z n, is proportional to 1 1 m ( z n ) and estimates of unknowns may be 

n 

constructed by maximizing some variant of this. If n(z) = 0 for a particular z then p(z) may not be 
estimated directly. 


5.3.6 Measurement Errors 

At various stages of the NASA work, derived quantities subject to measurement error are fed into further 
computations as central variates or explanatories. One can mention inclinations and altitudes in 
particular. Several types of procedures are available to handle this circumstance, for example grouping 
and likelihood methods (e.g.. Fuller [5.32]). NASA is presently handling the problem by grouping the 
values into cells or ignoring the uncertainty. We can expect a likelihood approach to be more efficient. 

In Section 5.3.3c above, the problem was set up as one of Poisson regression, as if the p kV were known. 
In fact they are subject to measurement error, having been estimated from the data as described in 
Section 5.3. 1 above. Schafer [5.33] provides an EM method for handling the presence of measurement 
error. Nakamura [5.34] develops a direct method of handling the Poisson regression case of interest. 
Measurement errors do have the advantage of providing a type of regularization to the estimates. 

Such extravariation and the estimation procedure employed need to be taken account of in uncertainty 
computations for the final estimates. 


5.3.7 Uncertainty Estimation 

A complete study will include estimates of the uncertainties of any estimates presented. One has to be 
particularly concerned about the uncertainties of the estimates of the size distribution and of the flux. 

The main methods available for this problem are maximum likelihood formulas, propagation of error 
(based on linear approximations), sample splitting, and jackknife and resampling methods (see Efron and 
Tibshirani [5.35] for a general survey). 
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Consider the equation (5.7) above. The function f(z) is random and p(z\s) is an estimate, i.e., two sources 
of variability feed into the estimate / . In the case of the estimate of Section c above and that the p k[l 
may be considered known, the standard errors of the /, are part of the standard output of a generalized 
linear model program. 

For the regularized case and when the estimates are determined via the EM method, Segal et al. [5.36] 
present a method for determining uncertainty if X may be considered known. Wahba [5.24] and 
O’Sullivan [5.19] also present methods for the case that p(z\s) is known. 

In the case that the sampling variation of p( 4) is appreciable the jackknife and resampling appear to be 
viable approaches to handle the two aforementioned sources of uncertainty. The bootstrap can be 
computer intensive so also studying the propagation of error, simple sample splitting, and jackknife 
approaches appears worthwhile. 


5.3.8 Goodness of Fit/Model Validation 

Contemporary statistics has a variety of tools to offer to the problems of goodness of fit and model 
validation: residual analysis, goodness of fit criteria, etc. Contemporary books on regression analysis, 
e.g., Cook and Weisberg [5.37], go into these in detail. 

NASA’s present approach is to compare cumulative flux rate and expected detection rate with a 
corresponding observed set of detections. (To quote a NASA handout, “Model parameters are then 
adjusted based on theory until its predictions conform to Haystack measurements and other data.”) The 
detections are assumed Poisson. This assumption seems reasonable since the detections rate is low. For 
some questions, one cannot ignore the presence of swarms or coherent groups of objects. A variety of 
point process models exist for such cluster processes. 

5.3.9 Forecasting 

The ORDEM96 program that NASA makes available requires the forecasting of some basic covariates, 
e.g. solar activity and launch rate. The implications of this, e.g. sensitivity to particular values, bears 
study. Section 5.3.2.C mentioned one method of analytically including covariates. Formal forecasting 
procedures exist that provide estimates of uncertainty in some cases (e.g. Harvey [5.38]). It appears 
sensible to proceed here by providing results for various scenarios and margins of safety. 
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6. Conclusions 

The Panel reached the following conclusions in regard to the issues posed: 
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Figure 6.1 Typical plots of coefficient of variation of the flux ( the standard error a of the flux divided 
by the flux) vs. total hours of observation for objects in different size intervals. To obtain the standard 
error of the flux, multiply the results on the y-axis by the mean flux over the time in terval. 

Issue 1 : “The number of observations relative to the estimated population of interest” 

The number of observations relative to the estimated population is a statistical question and the panel has 
taken a statistical approach to answering the question. An analytical procedure has been developed for 
calculating the number of observations required to estimate the target populations with prestated 
confidence bounds. The model takes into account known sources of error and assumes a family of 
distributions. 

Conditioned on a number of assumptions which are discussed in Appendix A, curves have been produced 
indicating the confidence interval for the population estimate as a function of the total number of 
observations available. Figure 6.1, based on a subset of the Haystack data collected in 1994 and 
representing 840 detections over a 97.22-hour observation period, shows the sort of behavior which 
might be expected as the number of observations increases. The curves in Fig. 6.1, parameterized as a 
function of debris particle size, show that confidence in the population estimate improves with increasing 
numbers of observations. These curves are illustrative of the type of results that can be obtained, but 
should not be taken as representing the actual values for the entire Haystack data set. 

The details of the confidence interval curves will change as the population upon which they are based 
changes. However, a statistical procedure, explained in Appendix A, is available for NASA to produce 
curves based on other data sets. 



Observation time (hr) 

Includes extra Poisson variation multiplier of 1.16 
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Issue 2 : “The inherent ambiguity of the measured dBsm and the inferred physical size” 

In the case of the conducting sphere, a given RCS in dBsm can be generated by as many as three different 
sphere diameters, producing ambiguity. The NASA size estimation model (SEM), although it follows the 
curve for average sphere RCS, was derived from many measurements at different aspect angles and 
frequencies on a sample set of 39 objects. The resulting translation of RCS to size is subject to 
uncertainty (not ambiguity), described by a distribution of size values about the mean, for a given RCS. 
Flux vs. size curves can be derived from Haystack data, along with estimates of their uncertainty, 
notwithstanding the variability of RCS with size, shape, and viewing angle. These distributions are taken 
into account in the development of the uncertainty bounds discussed under Issue #1 . 

Issue 3 : “The inherent aspect angle limitation in viewing each object and its relationship to the object’s 
geometry” 

A Haystack observation of an object is made at unknown aspects, representing samples from a possible 
47t steradians of angle about the axis of the object. The conversion of measured RCS to size, via the 
SEM curve, is based on the average overall aspect angles and the distribution about this average caused 
by object shape and aspect angle variation. The use of this distribution in the derivation of confidence 
limits takes into account the fact that RCS values from only a limited number of aspect angles (typically 
one) are observed in a Haystack measurement. 

Issue 4 : “Adequacy of the sample data set to characterize the debris populations potential geometry” 

The sample data set derived from a high-velocity impact experiment produced a distribution of shapes 
that is quite similar to the distributions of shapes obtained from other high-velocity impact and explosion 
experiments. Without direct physical sampling of the objects on orbit, this is the most useful evidence 
that the measured data set is representative of objects on orbit. Therefore the measured data set can 
properly be assumed adequate to characterize the space debris population’s potential geometry. 


7. Recommendations 

The panel makes the following recommendations to NASA: 

1 . A succinct summary document is needed which describes the debris problem, the Haystack 
measurement program, the procedure for converting RCS to size, and the resulting debris flux results 
obtained to date. This information is embedded in several existing reports, but needs to be integrated 
into a single document for use by those nonexperts in fields of radar and space debris. 

2. Continue to make use of modern statistical techniques in addressing the interpretation of the 
Haystack data. 
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APPENDIX A 


Some Analyses of Haystack RCS Data 

David R. Brillinger 
Statistics Department 
University of California, Berkeley 


A.l. Introduction 

A number of statistical methods for the analysis of radar cross section (RCS) data collected by Haystack 
are described and discussed in the body of this report. This appendix presents results of their practical 
application to some data collected by the facility. There are also simulation studies of some of the 
properties of the estimates. The work is preliminary based on a single data set and not implementing 
some of the features discussed in the report. 

The report by Barton et al [A.l] listed as a Priority 1 item “Calculate confidence interval in flux from 
Haystack data.” Obtaining such intervals as a function of object size, altitude, and orbital inclination is 
one goal of the analyses presented here and an example is provided in Figure A. 10. The problem appears 
to be usefully divided into that of estimating statistical densities or distributions and that of estimating a 
factor to multiply these up by to obtain overall counts and fluxes. The former quantities may be 
estimated in statistically stable fashion, while the latter factors are inherently variable. 

Various difficulties arise. In going from an object detected to an estimate of its size there are a number 
of sources of uncertainty. These include the measurement error of an RCS value (see Figure A.2a below) 
and the spread of the calibration distribution (see Figure A.4b below). This last is caused in part by the 
random angle of observation of the target’s orientation. Further, the data are a biased sample, the 
probability of detection of each particle is not the same. However, this probability may be estimated (see 
Figure A. 2b below). 

A basic problem may be set up formally as follows. Let p(s,z) refer to the joint distribution of the 
(standardized) size and measured RCS of an object selected from the population of interest. Let 
p s (5), p 7 (z) refer to the distribution’s marginals. The former is a quantity of particular interest, and 
observations from the latter are available. The conditional distribution p(z\s) may be estimated via a 
(statistical) calibration experiment and a study of the measurement error of an observed RCS. The 
distributions just defined are related via 

Pz (z) = j P(z\ s)p s ( s)ds (A. 1 ) 

In the case that there is an exact calibration relationship z = Ms), s =g(z), one has 

p(z\s) = 6(z~h(s )) 

where 8( ) denotes the Dirac delta function and (A.l) then gives 
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Pz(z) = p s {g{z))\g\z)\ or p s (s) = p z (h(s))\h\s)\ 


by standard changes of variables. 

An estimate based on these last will be referred to as a direct estimate, one simply goes from z to s values 
via s = g(z). Such estimates will be studied in the work below. So too will an indirect estimate based on 
the relationship (A. 1). Complications to implementing the procedures: The calibration relationship is 
not exact, because of the variable angle of viewing, nonequivalence of size, and RCS quantities; and the z 
values obtained via Haystack are subject to measurement error. 

Let us emphasize that a variety of assumptions are made in developing the procedures presented. These 
assumptions need to be checked to the extent possible. 

A.2. Some Haystack Data 

NASA provided the data studied; the results of a preliminary study are given in Matney [ A .2] . The data 
are for the year 1994, with a beam elevation between 74° and 76° and an azimuth between 85° and 95°. 
There were 840 detections over observation periods totaling 97.22 hours. For each case, one has RCS 
estimate, altitude estimate, inclination estimate, and time of observation amongst other quantities. 
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Figure A.l. Histogram of the 840 available RCS values. 
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Figure A.2. (a) Estimated standard error versus estimated RCS value, (b) Inverse of the estimated 

probability of detection for each object. 

Figure A. 1 provides a histogram of the unitless RCS values of the sample with the units decibels. It has a 
bell shape with some skewing to the right. The median is -12.82 dB. There are outlying values of 26.35 
and 36.96 dB. The distribution graphed is in part reflecting both the actual distribution of the debris 
sizes and the fact that Haystack has less sensitivity to objects of smaller sizes. 

There is unequal probability of detection, depending on how an object happens to pass through the beam 
and variations of the signal due to noise and scintillation of the target. Fig. A.2b graphs estimates of the 
inverse of this probability for the 840 cases. The factor is close to 1.0 for RCS above -10 dB 

The RCS values are subject to measurement error. The size of this is estimated as part of the detection 
process. Figure A.2a graphs the measurement standard error against the RCS value for the 840 cases. A 
smooth curve has been added for reference. The scatter of the errors is seen to increase generally once 
one gets below -15 dB. (In the simulation studies reported below it is simply taken to be 2.5 dB., the 
level to the right in Fig. A.2a.) 

Figure A.3a provides an estimate of the probability density function of the RCS values themselves. 

Figure A.3b presents a bias-corrected density estimate, weighting the observations inversely to their 
estimated probability of detection. The “correction” is noticeable below -10 dB and reflects the 
substantial variability below -25 dB. 
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Figure A3. (a) Estimated probability density of the RCS. 

(b) Bias-corrected estimate. 

A.3. The Calibration Experiment 

A laboratory study was carried out and a statistical calibration curve relating “known” size to measured 
RCS values obtained. This curve is monotonic and given by 

s-^\z! 71 if z>2.835 

5 = 6 ^4z/9ti 5 if z<.00122 < A - 2 ) 

s = g(z) in between 

where 5 and z are dimensionless size and RCS values respectively and the smooth function g(-) is given 
by its values at 23 points. The relationship is linear at the extremes (with different slopes), when 
expressed in logs. The relationship (A.2) will be written s = g(z) generally. 

The distribution of RCS for given size and random viewing angle has also been estimated. A specific 
formula for p(z\s) is reported in Bohannon et al. [A.3]. It is a mixture of a gamma and a lognormal. The 
sizes of the objects studied ranged from 1 cm to 30 cm in the largest dimension and they were illuminated 
with radiation of frequency 2.4 to 18 GHz. The original data were not available for the present analyses 
and their further study appears worthwhile. 
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Figure A.4. (a) Result of applying the calibration function (A.2) to the RCS values. 

(b) Results from simulating the distribution employing the density as estimated in 

Bohannon et al. [A.3], 

Figure A.4a provides the results of applying the calibration relationship (2) to the n = 840 observed RCS 
values, specifically the points (s- = g(Zj),Zj),j = 1 ,..,n are plotted. A noticeable feature is the bump in 
the transition region, around s = size/ X = .3. (These are dimensionless units. For Haystack A, = 3 cm.) 

Figure A.4b presents the results of simulating from the distribution p(z\s) as provided in Bohannon et al. 
[A.3], to obtain £. for given S- and then plotting the ( ) points. This is done to provide an 
indication of the variability in the relationship and it appears substantial. 

A.4. Direct Estimates 

Direct estimation comprises ignoring the observation bias and making use of the functional calibration 
relationship (A.2) rather than the distribution p(z\s). This type of estimate has been used in many studies 
and is usually referred to as the Size Estimation Model (SEM). 

Figure A.5a provides the estimated density of the estimated sizes, S- = g(z ] ) , computed from the 
observed distribution of RCS values via the relationship (A.2). There is a tail to the right with 
exponential looking falloff. One explanation for the falloff in intensity to the left is Haystack’ s 
insensitivity to objects in that region. 
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The square root is the variance stabilizing transformation for a density estimate and it is often convenient 
to work on this scale. Figure A.5b graphs the square root of the density estimate, together with 
± 2 standard error limits about it. The width of the interval is seen to be constant and the low density 
features have become better emphasized. 


Density estimate of directly estimated size 


(a) 


Density 


(b) 


sqrt 

(density) 



0.1 0.5 1.0 5.0 10.0 50.0 100.0 

size/lambda 


Sqrt(density estimate) of directly estimated size 



Figure A.5. (a) Estimated density of size simply using the relationship (A.2) on the 
observed RCS values, (b) Square root of the density estimate 
and ±2 standard error limits. 

In a variety of cases it is the cumulative number of objects above a given size that is of interest, rather 
than their density. Figure A. 6 graphs 

G(s) = Y J H(s J -s) ln (A.3) 

i 

the proportion of cases with estimated size greater or equal to the specified value, s. [Here H( ) is the 
Heaviside indicator function.] Expression (A.3) provides an estimate of G{s ) the probability that an 
object sampled has size greater or equal to s. 

The variance of G(s) may be estimated by G(s)[l — G(s)] / n and is indicated on Fig. A. 6 via 
± 2 standard error limits. (The bias is another matter.) The width of the limits is seen to increase as one 
moves to the larger sizes. The falloff of the curve is approximately linear from size about 0.3 up. The 
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linearity corresponds to an exponential distribution for log size above a cutoff point. Curvature present 
at the smaller sizes may be seen as resulting from the censoring of observations. 


One has to be concerned about another type of bias in the estimate (A.3). It occurs since . is but an 
estimate of the actual size S - . It is perturbed by the measurement error involved in Zj and by the spread 
of the distribution p(z\s). This issue may be addressed analytically, but is first studied via simulation. 

Direct estimate of proportion with size/lambda > s 



Figure A.6. Estimate of the proportion of values greater or equal to a specified size. 
Also provided are ±2 standard error limits. 


A.5. Some Simulation Studies 

As discussed in the main body of this report, one has to be concerned about the effects of the 
measurement error of the RCS values and about the uncertainty involved in the calibration operation. 
Figures A.7 and A.8 provide some information in this connection. 

To begin, the band Fig. A.7 presents is simply the ± 2 standard error bounds of Fig. A. 5b shaded in. 
This estimate was based on the values resulting from calibrating the observed RCS values using the 
relationship (A.2) and provides a basis for inferences concerning the debris sizes and risks of collision. 
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+-2 s.e. about sqrt(density estimate) & average of 25 simulations 



Figure A.7. The band provides ±2 standard error limits about the square root density 
estimate of Fig. A. 5. The line is the average of 25 simulations of sizes starting 
from the directly estimated sizes and adding both the RCS-size variability and 

measurement error. 


To study the direct estimate’s bias properties, these size values were taken as “truth,” and RCS values 
generated employing the distribution p(z\s) of Bohannon et al [A.3] together with an added measurement 
error (s.d. 2.5 dB). A density estimate is computed as before, i.e. based on newly estimated sizes. This 
whole procedure was repeated 25 times. The average of the 25 sqrt(density) estimates is the line of the 
figure. Notice the curve is outside the bounds at the lower sizes corresponding to a biased estimate in 
this region. This is the region with probability of detection not 1 and that phenomenon was not included 
in the simulations carried out. Generally speaking, the direct estimate appears plausible above, say, 

5 = 0 . 3 . 

To study the situation in further detail, observations were simulated from an exponential distribution, 
with decay parameter like that suggested by the present data. The simulated values were taken to 
represent the log sizes of actual debris. The band in Fig. A. 8a corresponds to +2 se limits around the 
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sqrt(density) estimate based on the exponentials for a sample of size of 840, and corresponds to the band 
of Fig. A.7. 
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Figure A.8. Results of simulation experiments taking log sizes from an exponential 
distribution. The band provides +2 standard error limits about the square root density 
estimate for exponen tial variate. The solid line is the density estimate having generated 
RCS values and further perturbed them by measurement error. 


Next, the 840 exponential values generated were put through the calibration distribution and 
measurement errors added as in the case of Fig. A.7. A density estimate was then computed based on 
these simulated values and corresponds to the line in Fig. A. 8a. The estimate stays reasonably in the 
band. Continuing the work, sample sizes of 2500 and 5000 were also considered; Figs. A. 8b and A. 8c 
provide the results. Now the effects of the calibration and measurement errors show themselves. The 
width of the band has gotten smaller and the estimate falls noticeably outside, particularly in the region 
of transition of the calibration curve (A.2). 
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In summary, for a sample size of 840, the direct estimate appears a useful estimate, but its bias shows 
itself in the case of larger sample sizes. 
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Figure A.9. Estimate of number of objects from the sample of 840 with size greater or equal to a 
specified size s. Also graphed are ±2 standard error bounds. 


A.6. Count and Flux Estimation 

The analyzed data set had a sample size of n = 840 detections made over a period of 97.22 hours. This 
corresponds to an overall detection rate of 8.64/hour. A longer period would have resulted in more 
detections. The total count is itself variable and this is one of the reasons that densities and standardized 
distributions were the quantities first studied. 

If G(s) denotes a cumulative distribution as above and n the number of detections, then the number of 
particles of size greater or equal to 5 that passed through the beam in that time period is nG{s) . 
Assuming Poisson variability, an estimated standard error is ^nG(s) . Fig. A.9 graphs the quantity 
nG(s) and +2 standard error limits for the data at hand. Section 7 will mention a correction for non 
Poisson variability. 
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In the case of future time intervals — as required for risk estimation, for example — the standard error of a 
predicted count will involve both the variability of the count anticipated and of the estimate G(s ) . 

Flux at altitude h = 900 km 



Figure A.10. Flux estimate for h = 900 km, Ah = 200 km, in units of count/m 2 -yr. 
Also included are ±2 standard error bounds. There were 512 objects in the sample. 


Consider next the problem of flux estimation. The data available are for a particular time period and the 
beam in a particular position. From this information we can make inferences concerning the total 
population of objects. Because of the spreading of the beam, flux will be estimated as a function of 
altitude. Supposing the beam width to be 6 , the staring angle to be Ot and the altitude range of interest to 
have width Ah , the area of the surface of the frustum of the beam will be approximately 

A A = InhAh tan.50 / sin a . (A.4) 

A flux estimate may now be obtained as 

n(h,s,T)/[AAT], (A.5) 

where n(h,s,T) is the number of objects, of size greater or equal s, detected in the altitude interval 
(h—.5AhJi+.5Ah) during the time interval of length T. For detections per unit time one would not 
divide by A A . 
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For the data being studied 0 = 0.05 and a = 75°. The uncertainty of the flux estimate will be that of 
n(h,s,T), i.e., an estimate of the standard error of the flux provided by ^n(h,s,T) / [AAT] under a 
Poisson assumption. 
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Figure A.ll. Spline-based representation of the density. Also included 
are ±2 standard error bounds and Fig. A.6. 


A.7. Maximum Likelihood Estimation 

Estimating unknown quantities via the method of maximum likelihood has been much studied in practice 
and is typically highly efficient. There are a number of ways to implement such an approach in the 
present situation. One follows. 

Suppose that the desired density function of the sizes is parameterized as 
p s (s ) = c(0) exp { £ 0 p B p (log 10 (^)) } 

P 
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where 0 = (0 1 ,...,0 iJ ) is to be estimated and the B p (s) denote polynomial spline functions. The 
multiplier c(0) makes the density integrate to 1 . The likelihood 


Pz(Zi)-P Z (z n ) 

may be expressed as a function of 0 via the relationship (A. 1) and 0 estimated by maximization and 
thence p s (s ) estimated. An estimate of P rob {.S' > ,v} may then be based on p s (s). 

A pilot study of the procedure has been carried out and the results may be seen in Figure A. 1 1 . The 
function p(z\s) was taken to be that of Bohannon et al [A. 3], i.e., the measurement error phenomenon was 
not included although it is clear how to do so. Spline functions with 5 nodes were employed. The 
± 2 standard error bounds were obtained via a logit transform and the jackknife procedure; see Efron 
and Tibshirani [A.4] for a discussion of the latter. The data are separated into 20 groups in the jackknife 
implementation. Compare Fig. A. 1 1 with the direct estimate of Fig. A. 6, and note a bowing at the 
smaller sizes in the new estimate, but remember that these computations are preliminary. 

The properties of such estimates need to be further studied, e.g., by simulation. Further, there are 
objective ways to estimate the number of nodes. One may proceed to flux estimates as above. 

The bias-corrected case has not yet been implemented, nor has the case with measurement error. 

A.8. Regularized Estimation 

Not yet implemented, but the formulae are available. 

A.9. Temporal Aspects 

In computing the uncertainty above, the overall count has been taken to be Poisson. This would be the 
case were the sequence of passage times a Poisson process. This assumption may be studied from the 
passage times detected. 

One has also to be concerned regarding the stationariness of the relationships over the time period of 
observation. Interpretation of the estimates is difficult when the environment is changing, and covariates 
need to be included in the model to handle that problem. 

Figure A. 12 presents information concerning the rate with which objects are passing through the beam as 
a function of time. The square roots of the estimated rates for the individual periods of observation are 
plotted. There are 94 periods. The square root is taken, as this often stabilizes the variance of count 
variates. Also graphed are ± 2 standard error limits and a line at a height corresponding to the overall 
rate. One does not notice a time trend; however more than Poisson variation appears present. The 
extravariation may be estimated by generalized linear model techniques (see McCullagh and Nelder 
[A.5]). The estimate of the dispersion parameter is 1 .340 for the counts of the various periods, in 
contrast to the value 1 for a Poisson. To incorporate the phenomenon, multiply up the Poisson standard 
errors by 1.16 here. 
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Figure A.12. Estimated rates for the various observation periods. The square 
roots of the rates are graphed as well as ±2 standard error limits. 

In Section 4 an exponential falloff of the size distribution was noted as the size increased. This can lead 
to useful simplifications and comparative studies. In working on problems of seismic risk assessment, it 
is often useful to consider so-called ^-values. Observed earthquake magnitudes are assumed to have an 
exponential distribution above some cutoff level. The lvalue corresponds to the decay rate of the 
exponential as the magnitude increases. Its size gives an indication of the relative numbers of small and 
large events. The temporal variation of b has been studied as an earthquake predictor. 

In the present case the /?- value would be estimated as 

b = l/(U-U 0 ) 

where U is the mean of the log s values above the cutoff U 0 . One sees that a small b means more large 
values, relatively. 

Figure A. 13 provides the ^-values for the individual observation periods. Notice an indication of lower 
values for some particular periods. Uncertainty needs to be included. 
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Figure A.13. b-values for the various observation periods taking a 

cutoff level of 0.2. 


A. 10. Summary and Discussion 

With the sample size at hand, n = 840, there appears no cause for immediate concern in employing the 
direct estimate, particularly since the maximum likelihood estimate itself is variable. With larger 
samples, the direct and indirect estimates may be expected to differ noticeably. What is occurring with 
the sample size of 840 is that, while the direct estimate is biased, it has smaller variability and so is 
effective overall. 

The computations were carried out via the statistical package Spins (see, e.g. Venables and Ripley [A.6]). 
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